#!/usr/bin/R

library("abc");

 set.seed( as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31) )


SIMULATIONS=10000
PROPORTION=0.10 ##1%... 
ACCEPTANCES=SIMULATIONS*PROPORTION;
print(ACCEPTANCES);
GENERATIONS=50


d<-read.table("emperical.txt")[,1]

scale=10

xbar=mean(d);
s=sd(d);
n=length(d);
u=d[1]; t=sqrt(var(d));  ## this is the ANALYTICAL data....
mean=u*1/t^2/(1/t^2+n/s^2)+xbar*n/s^2/(1/t^2+n/s^2);
sd=sqrt(((t^2*s^2/n)/(t^2+s^2/n))^2);
plotlims=seq(mean-scale*sd,mean+scale*sd,length.out=1000);




parameters=rep(mean(d),ACCEPTANCES)+rnorm(0,var(d),n=ACCEPTANCES);
last_parameters<-parameters;
last_weight<-rep(1/(ACCEPTANCES),ACCEPTANCES);
stat_obs=0;

for (i in 1:GENERATIONS) {

##here's the C code for weighting...
#float mvn_pdf(struct Rn *mu, struct svd *ed,struct Rn *true_mu

##                        mpi_progress.current_weights[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]=0;
#                        for (uint i=0;i<mpi_ranks*accepts_per_rank;i++)
 #                               mpi_progress.current_weights[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]+=mpi_progress.last_weights[i]*mvn_pdf(sample_mean,sample_ed,mpi_progress.mean->f[i]);
 #                               //if (mpi_progress.last_weights[i]*mvn_pdf(sample_mean,sample_ed,mpi_progress.mean->f[i])>1)
 #                       mpi_progress.current_weights[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]=mvn_pdf(sample_mean,ed,mean)/mpi_progress.current_weights[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]];


	##parameters<-parameters+rnorm(n=SIMULATIONS/10); #perturb the parameters with N(0,1)
	param_sim<-vector();
	true_sim<-vector();
	stat_sim<-vector();
	max_weight=max(last_weight);
	for (j in 1:SIMULATIONS) {
		while (rand_sample<-ceiling(ACCEPTANCES*runif(1))) { ##this picks which parameter to perturb...
			if (runif(1)<=last_weight[rand_sample]/max_weight) { 
				true_sim[j]<-parameters[rand_sample];
				param_sim[j]<-true_sim[j]+rnorm(0,2*var(last_parameters),n=1);##this is where we perturb... think this is where we double variance of last generation
				stat_sim[j]<-sum((d-rnorm(param_sim[j],var(d),n=20))^2) ##think this is where we hold it constant at the sample variance
				break;
			}
		}

	}
	rej<-abc(target=stat_obs,param=param_sim,sumstat=stat_sim,tol=PROPORTION,method="rejection");

	true_parameters<-true_sim[rej$region=="TRUE"];
	parameters<-param_sim[rej$region=="TRUE"] ##take the accepted parameters 

	##ONLY DO THIS ON THE ACCEPTED>>>
	current_weight<-rep(0,ACCEPTANCES);
	for (j in 1:(ACCEPTANCES)) {
		for (k in 1:ACCEPTANCES) {
	##print(length(current_weight))
	##print(k);
			current_weight[j]=current_weight[j]+last_weight[k]*dnorm(parameters[j],var(last_parameters),true_parameters[k]); 
		}
		current_weight[j]=dnorm(parameters[j],var(last_parameters),mean(d))/current_weight[j]; 
	}
	print(summary(rej));
	print(paste("Working on iteration ",i,sep=' '));	
	system(paste("rm -v hist",i,".png",sep=''));
	png(paste("hist",i,".png",sep=''));
	hist(parameters,xlim=c(6,8),ylim=c(0,18),freq=F,xlab=paste("K.S. n-1 vs n p-value: ",ks.test(parameters,last_parameters)[[2]],"\nK.S. n vs analytical p-value: ",ks.test(parameters,pnorm,mean,sd)[[2]],sep=''));
	lines(plotlims,dnorm(plotlims,mean=mean,sd=sd),col="red");
	print(ks.test(parameters,last_parameters));
	
	last_parameters<-parameters;

	##average distance of accepted parameters is: 
	print(mean(stat_sim[rej$region=="TRUE"]))
	dev.off();


last_weight<-current_weight;
}

system("rm -v hist.gif");
system(paste("convert -delay 50 -loop 0 hist{1..",GENERATIONS,"}.png hist.gif",sep=''));
